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Abstract. The numerical solution of high dimensional Vlasov equation is usually performed by 
particlc-in-ccU (PIC) methods. However, due to the well-known numerical noise, it is challenging to 
use PIC methods to get a precise description of the distribution function in phase space. To control 
the numerical error, we introduce an adaptive phase-space remapping which regularizes the particle 
distribution by periodically reconstructing the distribution function on a hierarchy of phase-space 
grids with high-order interpolations. The positivity of the distribution function can be preserved by 
using a local redistribution technique. The method has been successfully applied to a set of classical 
plasma problems in one dimension. In this paper, we present the algorithm for the two dimensional 
Vlasov-Poisson equations. An efficient Poisson solver with infinite domain boundary conditions is 
used. The parallel scalability of the algorithm on massively parallel computers will be discussed. 

Key words. Particle-in-cell (PIC) Methods, Adaptive Mesh Refinement, Phase-space Remap- 
ping, Numerical Noise, Vlasov-Poisson equation 

AMS subject classifications. 

1. Introduction. The Vlasov equation describes the dynamics of a species of 
charged particles under electromagnetic fields. In the electrostatic case, the normal- 
ized equation reads 

df 

-^+v-V^.f+ {-lYiE + E') ■ V„/ = 0, (1.1) 

where f{x,v,t) is the distribution function of the species in phase space {x,v) S 
jgid ^ jgid ^[^]^ d = 1,2,3. sisO for positive charges and is 1 otherwise. E and E'^ 
denote the self-consistent and the external electric field, respectively. This equation 
is the simplest model to study coUisionless plasmas and beam propagation which is 
of importance to controlled thermonuclear fusion and accelerator modeling. 

The Vlasov equation is a nonlinear hyperbolic equation in phase space so methods 
of solution can be guided by the well-established numerical analysis of classical partial 
differential equations. Accordingly, grid methods in fluid dynamics, such as transform 
methods, finite-volume methods, and semi-Lagrangian methods, can be employed. 
Operator splitting was successfully applied to the solution of the Vlasov equation by 
Cheng and Knorr |4j in 1970s. It reduces the solution of the multi-dimensional Vlasov 
equation to a set of one-dimensional advection problems, and therefore has become 
a widely used technique. Even with these well-established algorithms, performing 
high-dimensional simulations using grid methods is still a challenging task. The issue 
is the computational time and memory cost in dealing with the whole six dimensional 
phase space. With the advances of supercomputer, grids methods have achieved 
large development in the last decade. In semi-Lagrangian methods, Sonnendrucker 
et al. |3D] introduced the cubic spline method. Nakamura and Yabe [55] introduced 
the cubic interpolated propagation method. In finite-volume methods, Fijalkow |14) 
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presented the flux balance method. A similar idea is used in a high-order finite- volume 
method based on mapped coordinate by Colella, Dorr and Hittinger [5]. Filbet, 
Sonnendrucker, and Bertrand |16j proposed the positive and flux conservative scheme 
using the idea of limiter. 

A more widely used approach for the solution of the Vlasov equation is PIC meth- 
ods T^jT. In PIC methods, the particles, a Lagrangian discretization of the distribu- 
tion function, follow trajectories computed from the characteristic curves given by the 
Vlasov equation, whereas the self-consistent fields are calculated on a grid. Since the 
methods employ the fundamental equations without much approximation, it allows us 
to observe most of the physics in a plasma system with relatively few particles. How- 
ever, as with all other particle methods, PIC methods suffer from numerical noise 
such that they have difficulty in simulating some problems, e.g., the problem with 
large dynamic ranges in velocity space. To remedy this deficiency, there are usually 
two approaches. One method is the so-called 6f method [13l |28l [21], discretizing 
only the perturbation 5f with respect to an equilibrium state /o based on a particle 
method. The Sf method has been successfully used in realistic applications, e.g., mi- 
crotubulence in magnetic confined plasmas |22| . The limitation of this method is that 
it can only be applied to the problems which are close to equilibrium. An alternative 
approach is through periodically reconstructing the distribution function on a grid in 
phase space. Such remapping technique has been used in particle methods in fluid 
dynamics [5|f2], i.e., vortex methods and smoothed particle hydrodynamics (SPH), 
to maintain regularity of the particle distribution and thereby improve accuracy, but 
has much more limited use in PIC methods in plasma physics. It is worth mentioning 
that early work of Denavit 12J and more recent work of Vadlamani f3TJ and Yang 
[3j used the idea of remapping for PIC methods. However, They all used low order 
interpolation function which results in a first order method overall. 

We studied a high order remapping scheme to PIC methods for the solution of the 
one-dimensional Vlasov-Poisson equations early [32j . Meanwhile, we provided a local 
redistribution technique such that the positivity of the distribution function could be 
preserved after high-order remapping. The initial numerical experiments on a set of 
classical plasma problems in one dimension are very encouraging. Remapping sig- 
nificantly reduces the numerical noise and results in a more consistent second-order 
convergence rate in the electric field error. We also investigated the effects of integrat- 
ing mesh refinement to the uniform remapping. This is motivated by the observation 
that remapping, a numerical diffusive procedure, tends to create a large number of 
small-strength particles at the low density region of the distribution function. Mesh 
refinement has the potential to reduce this side effect. 

In this paper, we extend the algorithm to the solution of the two-dimensional 
Vlasov-Poisson equations. This includes the use of an efficient Poisson solver with in- 
finite domain boundary conditions for beam problems. High-dimensional simulations 
are very expensive with respect to memory usage. We perform the simulation on a 
parallel machine using domain decomposition in physical space. A scalable implemen- 
tation based on domain decomposition in phase space will be discussed. We consider 
two types of numerical tests: plasma problems including linear Landau damping and 
the two stream instability, and a beam problem based on the paraxial model |15) . 

The rest of the paper is organized as follows. In f|2j we first review the classical 
PIC methods for the Vlasov-Poisson equations. An efficient algorithm which solves 
the Poisson equation with infinite boundary conditions is described. Then we present 
the high-order and positive remapping on a hierarchy of locally-refined grids in two 
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dimensions. ^ discusses the parallel implementation of the algorithm. We show the 
numerical results in S|4j Conclusion and future research will be given at the end. 

2. Algorithms. 

2.1. PIC methods. PIC methods are based on the Lagrangian description of 
the Vlasov equation 

'li^^O, (2.1) 
dt ' ^ ' 

where the characteristics {X{t), V{t)) arc the solution of the equation of motion: 

^ = ^W' ^ = (-l)-^(^(^- 1) + E-'l^X, t)) (2.2) 

with initial conditions X{t = Q) = x and V^t = Q) ^ v. 

In the beginning, the distribution function is approximated by a collection of 
point particles, 

f{x, t>, < = 0) « ^ qk6{x - Xk)5{v - Vk), (2.3) 
fc 

where {xi^,Vk) is a initial particle location at the cell center of a grid in phase space 
(quite start), = f{xk,Vk,t = Q)hxhyh^^h^^ is the weight of a particle. Then each 
particle follows a trajectory described by the equation of motion, 

^=0, ^(t) = y,(t), '^J^it) = {-nEu{t)+El{t)), (2.4) 

where Xfe(t = 0) = a;^. and Vk{t = 0) v^- 

At any time that a smooth representation of the distribution function is required, 
we approximate the function with a collection of finite size particles, where the exact 
delta function is replaced by a smoothed delta function. That is, 

/(a;,-i;,i) «^gfcJ,Ja;-Xfc(t))(5,„(t;-Ffc(i)), i > 0. (2.5) 

fc 

The smoothed delta function satisfies 



JR2 



(2.6) 



and 

1 



where u is any interpolation function and e is the stencil size. Usually, the stencil 
size for the smoothed delta function in physical space is chosen as the same as 
the mesh spacing of the Poisson solver. The typical interpolation function for PIC 
methods is the first-order interpolation function 

.r(.)4'-'^' (2.8) 
I otherwise. 

The flow of a PIC scheme is 
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Assign particle charges on a grid in physical space, 

p{xj,t) =J2(lkSsAxj - Xk{t)), (2.9) 

where j G are the node index of the Cartesian grid in physical space. 
The grid size is chosen as the same as the stencil size of the smoothed delta 
function Sx- In the case of the first-order interpolation function, for each 
node j, the sum is restricted to the particles with \xj — Xk(t)\ < ex- 
Solve the Poisson equation on the grid with a second-order finite-difference 
method: 

- (A (j)). - - 2^ - (-1) pj + ^background (2. 10) 



d=0 



and 



V^J^e^^O^^ (2.11) 

'^background background charge density if applicable. With given 

boundary conditions, the discrete Poisson equation is usually solved by a 
fast Poisson solver, such as FFTs or multigrid methods. 
Interpolate the calculated field back to the particle locations with the same in- 



terpolation function in equation (2.7 1. It is worth mentioning that a different 
interpolation function will introduce self- force errors [7]. 
• Integrate the equation of motion numerically, for example, using the second- 
order Runge-Kutta method. 

2.2. Solving the Poisson equation with infinite domain boundary condi- 
tions. To model beam problems, the Poisson equation with infinite domain boundary 
conditions needs to be solved. We compute the solution using a new version of the 
James-Lackner method [2D] by McCorquodale et al. [53J [21] . This method solves two 
Dirichlet boundary problems plus a boundary to boundary convolution. 

We briefly describe the algorithm below. Assume is the support domain of the 
right-hand side p, we can solve the Poisson equation with infinite domain boundary 
conditions on a slightly larger domain Di > D with inhomogeneous Dirichlet bound- 
ary conditions. The boundary value can be calculated by Green's function convolution 
from the source p to the domain boundary dDi. The volume source p to boundary 
dDi convolution is relatively expensive, in particular for 31? problems. Instead of 
using a volume to boundary convolution, we can compute the boundary value by 
performing a boundary dDi to boundary dD2 convolution and solving another Pois- 
son equation on a domain D2 > Di > D with Dirichlet boundary condition. The 



procedure of James' algorithm is (Figure (2.1)) 



Step 1: Solve the Poisson equation on domain Di with homogeneous Dirichlet 
boundary conditions 

~A(t>i=p on Di, (t>i=0 on dDi. (2.12) 

Step 2: Calculate the surface charge on dDi 

dpi = ^ on dDi. (2.13) 
on 
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• Step 3: Perform a boundary to boundary convolution from dpi to dD2 

dct>2= I G{x-y)dpAy)dAy (2.14) 

JdDi 

using fast multiple methods [T7] . 

• Step 4: Solve another Poisson equation on domain D2 with inhomogeneous 
Dirichlet boundary conditions 

— A02 = p on D2, 02 = d(p2 on dD2- (2-15) 



du2 = G ■ dp 



U2 = du2 



wi = 



dp 



dui 

dn 



Aui = p 



"2 
AU2 = p 



Fig. 2.1. James Algorithm Left: Solve for ui on Di Middle: Calculate surface charge and 
perform convolution Right: Solve U2 on D2 



2.3. Particle Remapping. The convergence of particle methods for the one- 
dimensional Vlasov-Poisson equations has been investigated by Cottet and Raviart 
[2]. Their result shows that particle overlapping and regularization are important for 
the convergence of the methods. Specifically, the truncation error of a particle method 
is amplified by a time dependent exponential term. Based on Cottet and Raviart's 
work, we extend the error analysis to PIC methods [32] . Our result is one order higher 
in the truncation error. However, as in Cottet and Raviart's analysis, the truncation 
error is amplified by a time dependent exponential term. The analysis motivates 
the use of remapping technique, a widely used strategy in particle methods in fluid 
dynamics, to control the exponential error. The basic idea of remapping is simple. 
Since particles will gradually move away from the exact trajectories due to numerical 
error, we can reduce the displacement by periodically reproducing the distribution 
function f{x, v, t) on a grid by interpolation. A new set of particles, which are created 
from the grid representation, then replace the distorted particle distribution. The 
later step is identical to the initial step of PIC methods that we initialize the particle 
positions and weights in equation (2.3). The error due to remapping will depend on 
the order of the interpolation function. 

In the previous work [52], we successfully applied the remapped PIC method to 
the one-dimensional Vlasov-Poisson system. The remapping scheme was extended in 
three aspects compared with the standard scheme. First, we used high-order inter- 
polation functions which improve accuracy but do not preserve positivity. Second, 
we preserved the positivity of a high-order interpolation by redistributing the excess 
charge into its local neighborhood. The local redistribution algorithm is based on 
the mass redistribution idea of Chern and Colella [S] , which is first applied to enforce 
positivity preservation by Hilditch and Colella [18] . Third, instead of reinitializing 
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on a uniform grid, we reproduced the distribution function on a hierarchy of locaUy- 
refined grids. Remapping on a hierarchy of locally-refined grids significantly reduces 
the number of small-strength particles located at the tail of the distribution function. 
The high-order, positive, and adaptive remapping scheme in high dimensional phase 
space is described below. 

2.3.1. High-order Remapping. The overall accuracy introduced by remap- 
ping will be one order lower than the order of the interpolation function since we lose 
one order of accuracy in the evolution step. For example, the interpolation function 
with second-order accuracy only results in a first-order method overall. In this paper, 
we consider an interpolation function with third-order accuracy derived by Monaghan 
[25j . The function in one-dimensional can be expressed as 

< s = Jfi < 1 

l<s=^<2 (2.16) 
otherwise. 

The one-dimensional expression can be generalized to four dimensions by tensor prod- 
uct, 

3 

W4{xi-Xk)^Y[Wi{xi~x'i,h''), (2.17) 

where h'^ is the remapping mesh spacing in phase space, i and k denote the index for 
the cell-centered grid and the particles, respectively. 

This function, called a modified B-spline, conserves the total charge and represent 
a quadratic polynomial exactly. In addition, the first- and the second-order derivative 
of W4{x, h) are continuous. The smoothness property of this modified B-sphne is par- 
ticularly good for scattered data interpolation. However, as with all other high-order 
interpolation functions, W4 is not positivity preserving. An interpolation function 
without positivity might create nonphysical negative charge. This should be avoided 
in simulations. 

2.3.2. Positivity. The positivity preserving algorithm is based on the mass re- 
distribution idea of Chern and Colella [S] , first applied to enforce positivity preserva- 
tion by Hilditch and Colella [18]. In the algorithm, we redistribute the undershoot of 
cell i 

5/i = min(0,/r) (2.18) 
to its neighboring cells i + £ in proportion to their capacity ^ 

^i+, = max(0,/4,). (2.19) 
The distribution function is conserved, which fixes the constant of proportionality 

^-^+^+ ncJi^ (2-20) 

S Ci-l-fc 
fc^O 

for £ 7^ such that cell i + £ is a neighbor of cell i. Superscript n and n -I- 1 denote 
the interpolated value before and after redistribution, respectively. 

The drawback of this approach is that positivity is not guaranteed in a single 
pass. One might have to apply the method iteratively. In practice, however, we find 
a few iterations are sufficient. 
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2.3.3. Mesh Refinement. Mesh refinement is an attractive option in improv- 
ing the efficiency of phase-space remapping. The distribution function in phase space 
is inhomogeneous, for example, a MaxwelHan distribution in velocity space. When 
we represent the system by particles, it is best that we can have each particle car- 
ries a similar amount of weights. Remapping on a hierarchy of locally-refined grids 
is a good strategy for creating a such set of particles. From another point of view, 
remapping through interpolation is a numerically diffusive procedure. This results 
in a large number of small-strength particles near the tail of the distribution func- 
tion. The situation becomes worse as we apply remapping frequently. Remapping 
on a hierarchy of locally-refined grids, with a coarser grid covering the tail of the 
distribution function, can reduce the number of those small-strength particles. In the 
following, we present the algorithm of remapping with mesh refinement. In designing 
the algorithm, we have two guiding principles: the total charge should be conserved 
and the overall accuracy on the field needs to be maintained. 

Before explaining the algorithm, we introduce the definition of a composite grid. 
We define a hierarchy of cell-centered grids fig, where < £ < £max- ^i=o is the 
coarsest grid that covers the whole problem domain. The finer grids il.e>o are con- 



structed as a union of cell-centered rectangles (see Figure 2.2). The mesh spacing of 
each level is = hi_i/ri_i, where r^_i is the refinement ratio of level £—1. In four 
dimensions, he £ R'^ and e Z"'. The composite grid consists of valid grids at all 
levels, where a valid grid is defined as a region not overlain by a finer grid. That is, 

fie= £ (f^£\P/+l(f^£+l)), (2.21) 

e=o 

where P^^i is the operator projecting from level £ + 1 to level i. 

At the beginning, a set of particles are created from the cell center of the composite 
grid. In the remapping step, each particle first finds the valid cell in the composite 
grid it belongs to. One particle can only belong to a single valid cell. If the cell is 
far enough away from a coarse fine-interface, the charge can be interpolated on the 



grid as in equation (2.16). If the cell is near a coarse-fine interface such that the 
interpolation stencil intersects the coarse-fine interface, special care must be taken. 
First, we interpolate the charge on the surrounding cells as usual. After deposition, 
we know that not all deposited cells are valid cell. We need a further step to transfer 
the charge from invalid cell to valid cell. There are two cases depending on where the 
invalid cell is located. If the invalid cell is in a coarser level and it is covered by the 
valid cells of a finer level, we transfer the deposited charge from the coarser level to 
the finer level through interpolation. On another hand, if the invalid cell is outside the 
grid of the current level, the charge is transfered by projection. Figure [Z2| shows the 
algorithm in two dimensions. The four-dimensional case can be generalized easily. It is 
worth mentioning that we lose one order of accuracy in interpolating the coarser level 
charge into the finer level. However, since the coarse-fine interface is in co-dimension 
one, the expected accuracy in the field, e.g., second-order, will be preserved in Loo 
norm error. Our current implementation doesn't have time-dependent adaptivity. 
This feature can be incorporated by selecting some refinement criterion, for example 
each cell in phase space has similar number of particles. 

3. Parallel Implementation and Issues. The parallel implementation of the 
algorithm is straightforward based on domain decomposition in physical space. The 
physical space is decomposed into M patches. Given a parallel machine with N 



8 

















• 


• 


• 


• 






• 


• 


• 


• 






• 


• 


X 

• 


• 


• 






-< 

• 


• 








• 


• 


• 


• 


• 






• 


3 — 
• 





























































f-o- 


O 


-o-t 








— 1 

+: 


• 


• 


• 




• 




• 










— 1 — 


• 


• 


• 





































Fig. 2.2. Cross signs denote the particle locations, 
are denoted by filled circles and open circles, respectively, 
plots. Left: Particle is at the coarser level side. Right: 
denote the particle locations. 



The valid and the invalid deposited cells 
The refinement ratio is rg = (2, 2) in the 
Particle is at the finer level side. Cross 



processors, each patch is assigned to a processor cycUcally. Particles are assigned to 
patches according to their physical space positions. Using MPI, patches communicate 
with each other through ghost-cells and particles move between patches. This is the 
default implementation in Chombo software j27j . 

The current implementation is not a scalable algorithm (weak scaling) because of 
decomposition in physical space only. The potential issue is that when the problem size 
increases, since the number of processors is scaled in proportion to the problem size in 
physical space, the computational time and memory usage will increase in proportion 
to the problem size in velocity space. In the worst case, the processor will be out of 
memory. An alternative implementation is based on domain decomposition in phase 
space. In this implementation, it might be the case that particles belong to the same 



cell in physical space (equation (2.9)) are distributed on different processors. Since the 
Laplacian operator is linear, we can choose to solve the Poisson equation separately 
on different processors. The total fields are then obtained by MPLAUreduce. 

4. Numerical Tests. We demonstrate PIC methods with adaptive phase-space 
remapping on a set of classical plasma and beam problems in two dimensions, in- 
cluding linear Landau damping, the two stream instability, and beam propaga- 
tion in the paraxial model. For satisfying the overlapping condition, we choose 

£x/hx = £y/hy = 2, where = {£x-,£y)- 

We use Richardson extrapolation for error estimate. If E is the electric field 
computed with the initial phase space discretization (/i^,, hy^ h.^^, ^vy) ^^'^ integration 

step size At, and E computed with (2/ia,, 2/ij,, 2/i„^, 2/i^^) and 2At, the relative 
solution error in direction d is defined as 

eS-l^d-^fl- (4.1) 
q is the order of the method and is calculated by 

q = minlog2 ff^y (4.2) 
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4.1. Linear Landau Damping. 

damping is 

1 



Tlie initial distribution for linear Landau 



foix,y,v^,Vy) = — exp(-(u^ + v )/2){l + acos{kxx) cos{kyy)), 
Zn 



(4.3) 



0.5 and v„ 



6.0. The physical domain is {x,y) G 



where a = 0.05, = k 
[0,L = 2Tr/kx] X [Q,L = 2tt /ky] with periodic boundary conditions. Particle charges 
with strength less than 1.0 x 10^^ are ignored. In the simulation, we apply remapping 
every 5 PIC time steps. 

In the first test, we are interested in the evolution of the amplitude of the electric 
field. According to Landau's theory, the electric field is expected to decrease expo- 
nentially with damping rate 7 = —0.394. The behavior of exponential decay has been 
observed by many other authors, mostly calculated by grid methods [2^ [T51 [TUl [TT] . 

We initialize the problem on two levels of grids with base level at hx = hy = 
L/32, hy^ = hy = Wmax/16. The velocity space is refined on sub-domain v S [—3, 3] x 
[—3,3] with a refinement ratio 2. The PIC step size is dt ~ 1/8. We compare the 
simulation with and without remapping in Figure (4.1 ). In the case with remapping, 
the computed damping rate is very close to the theoretical value. The simulation 
without remapping fails to track the exponential decay. 





(a) with remapping 



(b) without remapping 



Fig. 4.1. The amplitude of the electric field for 2D linear Landau damping problem. Scales 
{hx, hv^ ) above denote the particle grid mesh spacing at the base level, where hx = hy and hy^ = h^^^ . 
With remapping, the computed damping rate is very close to the theoretical value 7 = —0.394. 
Without remapping, the simulation fails to track the exponential decay after a few damping circles. 



In the second test, we compare the electric field errors and corresponding con- 



vergence rates with and without remapping in Figure (4.2) and (4.3 1. We see that 



remapping significantly reduces the electric field errors and improves their correspond- 
ing convergence rates. 

4.2. The Two Stream Instability. The initial distribution for the two stream 
instability is 



fo{x,y,Vx,Vy) 



1 

12^ 



exp(-(w^ + vl)/2){l + a cos{kxx)){l + 5vl), 



(4.4) 



where a = 0.05, k^ — 0.5 and Wmax = 9-0. The physical domain is {x, y) e [0, L = 
2'K/kx\ X [0, L = 211 /ky] with periodic boundary conditions. Particle charges with 
strength less than 1.0 x 10^^ are ignored. As in linear Landau damping problem, we 
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(a) Loo errors of Ex 



h,;ui28,hVh:>4 




(b) Convergence rates 



Fig. 4.2. Error and convergence rate plots for 2D linear Landau damping problem with 
remapping. Scales {hx,hy^) above denote the particle grid mesh spacing at the base level, where 
hx = hy and hy^ = hy^. Second-order convergence rates are obtained, (a) the Loo norm of the 
electric field errors on three different resolutions, (b) the convergence rates for the errors on plot 
(a). 
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' h^=L!/64, li„ =h ' „/3^ 
hx=LV128, h„^=h^3,/64 




(b) Convergence rates 



Fig. 4.3. Error and convergence rate plots for 2D linear Landau damping problem without 
remapping. Scales {hx,hy) above denote the particle grid mesh spacing at the base level, where 
hx = hy and h^^ = hy^ . The errors without remapping are much larger than the case with remapping 
(see Figure J^.gp j. (a) the Lao norm of the electric field errors on three different resolutions, (b) 
the convergence rate for the errors on plot (a). 



apply remapping every 5 PIC time steps. The velocity space is refined on sub-domain 
V E [—4.5,4.5] X [—4.5,4.5] with a refinement ratio 2. 

We compare the electric field errors and their convergence rates with and without 
remapping as before. Figure ( 4.5 ) shows the Loo norm of the errors at the case without 
remapping in three different resolutions. The corresponding convergence rates are 
shown on the right of the error plots. Second-order convergence rates are lost at the 
early time of the simulation. Comparing with the results with remapping in Figure 
(4.4), we see that remapping extends the second-order convergence rates to longer 
times. 

We also compare the projected distribution function on plane (x^v^) at the same 
instant time t = 20 by both methods in Figure (4.6 ). For visualization purpose, in the 
case without remapping, we interpolate the particle-based distribution function to a 
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grid in phase space. We see that the classical PIC method results in a noisy solution 



see Figure (4.6(b))). Figure (4.6(a) I shows the distribution function computed by 



the PIC method with remapping. 




h,=L/64, h„=h„„J32 
h,=U128, h„=h^,,/64 




(a) Loo errors of 



(b) Convergence rates of the errors on the left 



Fig. 4.4. Error and convergence rate plots for the two stream instability with remapping. 
Scales {hx,hy) above denote the particle grid mesh spacing at the base level, where h^ = hy and 
hv^ = hvy ■ The PIC step size is dt = 1/8 at the lowest resolution, (a) the Lex, norm of the electric 
field errors on three different resolutions, (b) the convergence rate for the errors on plot (a). 
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(a) Loo errors of Ex 



(b) Convergence rates of the errors on the left 



Fig. 4.5. Error and convergence rate plots for the two stream instability without remapping. 
Scales (hxyhv) above denote the particle grid mesh spacing at the base level, where hx = hy and 
hv^ = hy^. The PIC step size is dt = 1/8 at the lowest resolution, (a) the Loo norm of the electric 
field errors on three different resolutions, (b) the convergence rate for the errors on plot (a). 



4.3. Semi-Gaussian Beam. The paraxial model is an approximation to the 
steady-state Vlasov-Maxwell equation in three dimensions. The K-V distribution is 
a measure solution of the paraxial model. Given an arbitrary initial distribution, we 
can focus a beam with the same matching forces for the K-V beam using the concept 
of equivalent beam. Here, we consider an initial semi-Gaussian beam focused by an 
uniform electric field using the concept of equivalent beam. The model has been 
considered by many authors [221 UHl ITT] . 

In the test, the beam is composed of ionized potassium. The physical parameters 
are the following: current / = beam velocity Vb — 0.63 x lO^m/s, and the radius 
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(a) with remapping at time t = 20 



(b) without remapping at time t = 20 



Fig. 4.6. Comparison of F{x,Vx), the distribution function projected on space (x,Vx), at 
the same instant of time t = 20 with (Left) and without remapping (Right) for the two stream 
instability problem. The projected value is F{x,Vx) = f{x,y,Vx,'Vy)dydVy. The grid-based 

distribution function is obtained by reproducing the particle-based distribution function through a 
second-order interpolation. We initialize the distribution function on two levels of grids, with base 



on sub-domain v £ [—4.5, 4.5] X [—4.5,4.5]. The classical PIC method results in a noisy solution with 
large errors in maximum. Both numerical noise and errors in maximum are significantly reduced by 
using remapping. The negative minimum in the case without remapping is a superficial effect due 
to project 4D data to 2D using a high-order interpolation for visualization purpose. 



of the beam a = 0.02to. We choose the tune depression rj = 1/2. For the normaUzation 
of the paraxial model, we refer to the work of Filbet and Sonnendrucker 15 . We use 
normalization parameters (xq^vq) — (a, ^f^i^)- This results in the normalized Poisson 
system 

-A(j)^ [ fdv, -\/(j)^E, (4.5) 

with initial semi-Gaussian distribution 

,t ^ / ^^^exp(-(i;2+«2)/2), + < 1 

I 0, otherwise 



and the external matching field 



E''{x,y,t) = - — {xe^ + yey). 
77^ 



(4.7) 



For the numerical parameters, we choose {hx — \Lx\/128, hy — \Ly 1/128) and {h^^ = 
w/128,/i,^ = Wmax/128), where {L,,Ly) = (-10,10) and w = 10. The PIC 
time step is dt = 0.00052925. 



Figure (4.7 1 shows the projection of the distribution function on planes {x,Vx) 
with and without remapping, respectively. The simulation with remapping gives a 
well-resolved result which preserves the positivity of the distribution function. IVIean- 
while, we show the root mean square (RIVIS) quantities of the semi-Gaussian beam 
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in Figure (4.8). Although the RMS quantities of the semi-Gaussian beam are os- 
cillatory, they remain close to the quantities of the associated K-V beam and they 
converge as the resolution increases. As mentioned by the other authors |10| . the 
oscillatory behavior is due to the fact that the semi-Gaussian is not exactly a steady 
state distribution. 



5. Conclusion. In this paper, we have presented the adaptive remapped PIC 
method to the high-dimensional Vlasov equation and demonstrated in linear Landau 
damping, the two stream instability, and the beam propagation problems. The new 
method reduces the numerical noise significantly. There are two extensions of the 
current research. The first will be the development of a scalable algorithm based 
on domain decomposition in phase space. The second will be the introduction of 
time adaptivity to the current algorithm that the hierarchy of locally-refined grids are 
dynamically created from the particle distribution at every remapping step. 
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PIC method results in a noisy solution with large errors in maximum. 
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